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Abstract. We propose a method to study lattice QCD at finite temperature 
(T) and chemical potential (/z). We test the method and compare it with the 
Glasgow method using n/=4 staggered QCD with imaginary /i. The critical 
endpoint (E) of QCD on the Re(//)-T plane is located. We use n./=2+l dynam- 
. ical staggered quarks with semi-realistic masses on Lt=4 lattices. Our results 

Tj- ' are based on O(I0 3 — I0 4 ) configurations. 
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1. Introduction 

QCD at finite T and/or fi is of fundamental importance, since it describes relevant 
features of particle physics in the early universe, in neutron stars and in heavy ion 
collisions. Extensive experimental work has been done with heavy ion collisions 
at CERN and Brookhaven to explore the pb-T phase boundary (cf. [0). Note, 
that past, present and future heavy ion experiments with always higher and higher 
energies produce states closer and closer to the T axis of the fi-T diagram. It 
is a long-standing question, whether a critical point exists on the fi-T plane, and 
particularly how to tell its location theoretically [ |^| . 

Let us discuss the /i=0 case first. Universal arguments [ [| and lattice sim- 
ulations [ [l| indicate that in a hypothetical QCD with a strange (s) quark mass 
(m s ) as small as the up (u) and down (d) quark masses (m u ,d) there would be a 
first order finite T phase transition. The other extreme case (n/=2) with light u/d 
quarks but with an infinitely large m s there would be no phase transition, only 
a crossover. Observables change rapidly during a crossover, but no singularities 
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appear. Between the two extremes there is a critical strange mass {ml) at which 
one has a second order finite T phase transition. Staggered lattice results on L t =4 
lattices with two light quarks and m s around the transition T (nf =2+1) indicated 
[ H that m c s is about half of the physical m s . Thus, in the real world we probably 
have a crossover. 

Returning to a non- vanishing /i, one realizes that arguments based on a variety 
of models (see e.g. [ ||, [?], ||) predict a first order finite T phase transition at large 
fi. Combining the /i = and large [i informations an interesting picture emerges 
on the \x-T plane. For the physical m s the hrst order phase transitions at large 
fi should be connected with the crossover on the /i = axis. This suggests that 
the phase diagram features a critical endpoint E (with chemical potential [Ie & n d 
temperature Is), at which the line of hrst order phase transitions (/i > [Ie and 
T < Te) ends [ ||. At this point the phase transition is of second order and long 
wavelength fluctuations appear, which results in (see e.g. [ ||) consequences, similar 
to critical opalescence. Passing close enough to (he,Te) one expects simultaneous 
appearance of signatures which exhibit nonmonotonic dependence on the control 
parameters [ [|, since one can miss the critical point on either sides. 

The location of E is an unambiguous, non-perturbative prediction of QCD. No 
ab initio, lattice analysis based on QCD was done to locate the endpoint. Crude 
models with m s = oo were used (e.g. [ ^]) suggesting that \ie ~ 700 MeV, which 
should be smaller for finite m s . The goal of our exploratory work is to propose a 
new method to study lattice QCD at finite /i and apply it to locate the endpoint. 
We use full QCD with dynamical nf =2+1 staggered quarks. 

QCD at finite /i can be given on the lattice [ [H)|]; however, standard Monte- 
Carlo techniques fail. At Re(/i)^0 the determinant of the Euclidean Dirac operator 
is complex, which spoils any importance sampling method. 

Several suggestions were studied in detail to solve the problem. We list a few 
of them (for a recent review see Ref. [ [ll| ) . 

In the large gauge coupling limit a monomer-dimer algorithm was used [ 
For small gauge coupling an attractive approach is the "Glasgow method" [ 
in which the partition function is expanded in powers of exp(/i/T) by using an 
ensemble of configurations weighted by the /Lt=0 action. After collecting more than 
20 million configurations only unphysical results were obtained: a premature onset 
transition. The reason is that the fi=0 ensemble does not overlap sufficiently with 
the states of interest. Another possibility is to separate the absolute value and the 
phase of the fermionic determinant and use the former to generate configurations 
and the latter in observables [ |l4| . The curvature of the fi — T phase diagram at 
/i = was determined by a stochastic calculation of the derivatives of the fcrmion 
determinant in [ . 

At imaginary jj, the measure remains positive and standard Monte Carlo tech- 
niques apply. The canonical partition function can be obtained by a Fourier trans- 
form [ [H| 0. In this technique the dominant source of errors is the Fourier 
transform rather than the poor overlap. One can also use the fact that the par- 
tition function away from the transition line should be an analytic function of /i, 
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and the fit for imaginary /i values could be analytically continued to real values of 
fi [ ^) . At T sufficiently above the transition, both real and imaginary /i can be 
studied by dimensionally reducing QCD 
help studying the problem [ M . 
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Hamiltonian formulation may also 




Fig. 1. Schematic difference between the present and the Glasgow methods. 

An attractive approach to alleviate the problem is the "Glasgow method" (see 
e.g. Ref. [ ^3|) in which the partition function (Z) is expanded in powers of 
exp(/i/T) by using an ensemble of configurations weighted by the /i=0 action. After 
collecting more than 20 million configurations only unphysical results were obtained. 
The reason is that the [i=0 ensemble does not overlap enough with the finite density 
states of interest [ ^l) . 

We propose a method to reduce the overlap problem and determine the phase 
diagram in the /i-T plane (for details see [ |22|). The idea is to produce an en- 
semble of QCD configurations at /i=0 and at the transition temperature T c . Then 
we determine the Boltzmann weights [ Q of these configurations at (i ^ and 
at T lowered to the transition temperatures at this non-vanishing fi. Since tran- 
sition configurations are reweighted to transition ones a much better overlap can 
be observed than by reweighting pure hadronic configurations to transition ones [ 
[l3fl . Since the original ensemble is collected at [i=0 we do not expect to be able to 
decsribe the physics of the large fx region with e.g. exotic colour superconductivity. 
Fortunately, the typical /i values at present heavy ion accelerators are smaller than 
the region we cover (for the freeze-out condition see e.g. [ psf ). 
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After illustrating the applicability of the method we locate the critical point of 
QCD. (Multi-dimensional reweighting was successful for determining the endpoint 
of the hot electroweak plasma [ e.g. on 4D lattices.) 

2. Overlap improving multi-parameter reweighting 

Let us study a generic system of fermions ip and bosons 0, where the fermion 
Lagrange density is , tjjM((j))ip. Integrating over the Grassmann fields we get: 



where a denotes a set of parameters of the Lagrangian. In the case of staggered QCD 
a consists of (3, m q and \x. For some choice of the parameters a—a a importance 
sampling can be done (e.g. for Re(/z)=0). Rewriting eq. (II]) 



We treat the curly bracket as an observable (measured on each configuration) and 
the rest as the measure. Changing only one parameter of the ensemble generated at 
ao provides an accurate value for some observables only for high statistics. This is 
ensured by rare fluctuations as the mismatched measure occasionally sampled the 
regions where the integrand is large. This is the overlap problem. Having several 
parameters the set «o can be adjusted to get a better overlap than obtained by 
varying only one parameter. 

The basic idea of the method as applied to dynamical QCD can be summarized 
as follows. We study the system at Re(/i)=0 around its transition point. Using a 
Glasgow-type technique we calculate the determinants for each configuration for a 
set of /i, which, similarly to the Ferrenberg-Swendsen method [ pl| . can be used for 
reweighting. The average plaquette values can be used to perform an additional 
reweighting in (3. Since transition configurations were reweighted to transition ones 
a much better overlap can be observed than by reweighting pure hadronic configu- 
rations to transition ones as done by the Glasgow-type techniques. The differences 
between the two methods are shown in Figure 0. Moving along the transition line 



We have directly tested these ideas in nf—A QCD with m q =0.05 dynamical 
staggered quarks. We first collected 1200 independent V=4-6 3 configurations at 
Re(/i)=Im(yLt)=0 and some /3 values and used the Glasgow-reweighting and also our 
technique to study Re(/x)=0, lm(/i)^0. At Re(/i)=0, lm(/i)^0 direct simulations 
are possible. After performing these direct simulations as well, a clear comparison 
can be done. Figure || shows the predictions of the three methods for the average 
quark condensates at /3=5.085 as a function of Im(^). The predictions of our method 
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was also suggested by Ref . [ |l7j . 
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Fig. 2. tj)ip as a function of Im(/x), for direct results (squares), our technique 
(crosses) and Glasgow- type reweighting (dots). 



agree with the direct results, whereas for larger Im(/i) the predictions of the Glasgow 
method are by several standard deviations off. We expect that our method can be 
applied at Re(/i)^0. 

3. The endpoint of n f = 2 + 1 QCD 

In QCD with n/ staggered quarks one should change the determinants to their J1//4 
power in our two equations. Importance sampling works also in this case at some /3 
and at Re(/i)=0. Since detM is complex an additional problem arises, one should 
choose among the possible Riemann-sheets of the fractional power in eq. (J2J) . This 
can be done by using [ the fact that at H=fi w the ratio of the determinants is 1 
and it should be a continuous function of /x. 

In the following we keep fi real and look for the zeros of Z for complex (3. At 
a first order phase transition the free energy oc logZ(/3) is non-analytic. A phase 
transition appears only in the oo limit, but not in a finite V. Nevertheless, Z 
has zeros at finite V, generating the non-analyticity of the free energy, the Lee- Yang 
zeros [|2(|. These are at complex parameters (e.g. (3). For a system with first order 
transition these zeros approach the real axis as V— ► oo by a 1/V scaling. This 
V— > oo limit generates the non-analyticity of the free energy. For a system with 
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Fig. 3. Im(/3o°) as a function of the chemical potential. 



crossover Z is analytic, and the zeros do not approach the real axis as V— *• oo. 

At T^O we used L t =A, i s =4,6,8 lattices. T=0 runs were done on 10 3 - 16 
lattices. m Ui( j=0.025 and to s =0.2 were our bare quark masses. At T ^ we 
determined the complex valued Lee- Yang zeros, (3$, for different V-s as a function 
of fjb. Their oo limit was given by a 0o(V) — f3fi° + Q/V extrapolation. We used 
14000, 3600 and 840 configurations on L s =4fi and 8 lattices, respectively. Im(/3g°) 
is shown on Figure |as a function of /i. For small fi-s the extrapolated Im(/?Q°) is 
inconsistent with a vanishing value, and predicts a crossover. Increasing /i the value 
of Im(/3o°) decreases, thus the transition becomes consistent with a first order phase 
transition (overshooting is a finite V effect). Our primary result is fiend = 0.375(20). 

To set the physical scale we used a weighted average of Ro, m p and -/ex. Note, 
that (including systematics due to finite V) we have (i?o • rn v ) — 0.73(6), which is 
at least twice, m u> d is at least four times as large as the physical values. 

Figure ^ shows the phase diagram in physical units, thus T as a function of /i^, 
the baryonic chemical potential (which is three times larger then the quark chemical 
potential). The endpoint is at T E = 160 ± 3.5 MeV, (i E = 725 ± 35 MeV. At ^=0 
we obtained T c = 172 ± 3 MeV. 
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Fig. 4. The T-/i diagram. Direct results are given with errorbars. Dotted line shows 
the crossover, solid line the first order transition. The box gives the uncertainties 
of the endpoint. 



4. Conclusion 

We proposed a method -an overlap improving multi-parameter reweighting technique- 
to numerically study non-zero fi and determine the phase diagram in the T-fi plane. 
Our method is applicable to any number of Wilson or staggered quarks. As a direct 
test we showed that for lm(/x)^0 the predictions of our method are in complete 
agreement with the direct simulations, whereas the Glasgow method suffers from 
the well-known overlap problem. We studied the fi-T phase diagram of QCD with 
dynamical n/=2+l quarks. Using our method we obtained Tb=160±3.5 MeV and 
/i£=725±35 MeV for the endpoint. Though /ig is too large to be studied at RHIC 
or LHC, the endpoint would probably move closer to the [1=0 axis when the quark 
masses get reduced. At ^=0 we obtained T c =172±3 MeV. More work is needed to 
get the final values by extrapolating in the R-algorithm and to the thermodynamic, 
chiral and continuum limits. The details of the presented results can be found in [ 
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